Construction of accurate Kohn-Sham potentials for the lowest states of 
the helium atom: Accurate test of the ionization-potential theorem 
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Accurate local Kohn-Sham potentials have been constructed for the ground Is 2 1 S state and, in 
particular, for the lowest triplet ls2s S state of the helium atom, using electron densities from 
many-body calculations and the procedure of van Leeuwen and Baerends (Phys. Rev. A49, 2138 
(1994)). The resulting Kohn-Sham orbitals reproduce the many-body densities very accurately, and 
furthermore we have demonstrated that the negative of the energy eigenvalue of the outermost elec- 
tron orbital agrees with the corresponding ionization energy with extreme accuracy. The procedure 
is also applied to the Hartree-Fock density of the ls2s 3 S state, and the Kohn-Sham eigenvalue of 
the 2s orbital is found to agree very well with the corresponding Hartree-Fock eigenvalue, which is 
the negative of the ionization energy in this model due to Koopmans' theorem. The results for the 
ls2s 3 S' state clearly demonstrate that there is no conflict between the locality of the Kohn-Sham 
potential and the exclusion principle, as claimed by Nesbet (Phys. Rev. A58, R12 (1998)). 

PACS numbers: 31.15Ew, 31.15Pf, 02.30Sa 



I. THE KOHN-SHAM MODEL 

According to the Hohenberg-Kohn (HK) theorem 0, , the energy of any electronic system 
can be expressed as a functional of the electron density, p(r), 

E[p] = FhkW + / drp(r)v{r), (1) 

where v(r) is the external potential and -Fhk[/>] is the universal HK functional, which in the 
constrained-search formulation is 0, |j| 

Fbk\p] =min*_ v (*|f + W\*). (2) 

Here, T is the kinetic-energy and W the electron-electron-interaction operators of the system (in 
atomic units), 



N N 

T = - 

i=l i<j 



N N 



The wave function, 'J, is normalized and belongs to the Sobolev space and the 

corresponding functional is defined for all A^-representable densities Q. The ground-state energy 
of the system is obtained by minimizing the energy functional over these densities 

E = mm p ^ N E[p] = E[p Q }. (4) 

This leads to the Euler-Lagrange equation 

— — — h v(r) = p, (5) 
5p(r) 

where p is the Lagrange parameter for the normalization constraint, J dr p(r) = N. 



* ingvar.lindgren@fy.chalmers.se 
t f3asos@fy.chalmers.se 

* gu99frmo@dd.chalmers.se 



2 



In the Kohn-Sham (KS) model the interacting system is replaced by a system of noninteracting 
electrons, moving in the local KS potential, t | Ks( r *) |2]i 

[-l v2 +%s«]^,W=^,(r). (6) 

The energy functional for this system is 

E KS [p] - T KS [p] + J drp(r)v KS (r), (7) 

where 

JV 

p(r)=£>(r)| a (8) 
i=i 

is the electron density. The kinetic-energy functional is 

T KS [p] =mm»^ p ($\f\9), (9) 

where <!> = det{0i} is a single Slater-determinantal wave function. Minimizing this functional leads 
to the Euler-Lagrange equation 

-j^ + v^ s {r)= l x. (10) 
Comparing with Eq. , leads - apart from an additive constant - to the relation 

VKs{r) = ^pW~^pW + v(r) - 

The Hohenberg-Kohn-Sham model was originally proven for the ground state but it was demon- 
strated by Gunnarsson and Lundqvist Q that it is valid also for the lowest state of a given 
symmetry. Later it has been shown to hold also for more general excited states |9j, [lOj . 

Although the form of the KS potential is generally not known, it can be constructed with 
arbitrary accuracy in cases where the electron density is known from other sources, e.g., from 
experiments or from ab initio calculations. Essentially two schemes have been developed for this 
purpose, by Zhao and Parr [TTllT2l[T^| and by van Leeuwen and Baerends [bij. respectively. 

The KS orbitals were originally assumed to have no other physical significance than generating 
the exact electron density, but it was later found by Perdew et al. 0, 0, and independently 
by Almbladh and Pendroza 

El 

that the eigenvalue of the outermost electron (with opposite sign) 
equals the ionization energy of the system. Perdew et al. have shown that considering densities 
that integrate to non-integrals, 



M 



J drp(r), (12) 



the theorem holds in the range N — 1 < M < N and hence when this number approaches N from 
below. This condition is known as the ionization-potential theorem |17| . 

The validity of the ionization-potential theorem has been challenged by Kleinman fl9j with 
counterarguments supplied by Perdew and Levy |17|- A number of numerical verifications of the 
theorem have been performed in the past [II Ifllll El El El El El E3 , generally with low 
or moderate accuracy due to problems in representing the density accurately using an analytical 
basis set [ilj. |23|. In the present work we use a numerical basis set, which has made it possible to 
demonstrate the validity of the theorem with much higher accuracy than in any previous calculation 
known to us |45j . 

The construction of the KS potential from the electron density has so far mainly been performed 
for atomic and molecular ground states. Recently, however, Harbola [27| has constructed the 
potential for the first excited singlet state, ls2s 1 S, of the helium atom. 
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Our primary goal for the present work has been to construct the KS potential for an excited 
triplet state, which, as far as we know, has not been done before. It has been rigorously shown 
that the KS potential is under general conditions strictly local 0, OJ I3, |2jj , and this has also 
been demonstrated in some of our previous works [sS, E0, E3] . Nevertheless, this fact has been 
disputed in several papers by Nesbet |33ll33 . l35| . who claims that the locality condition is in conflict 
with the exclusion principle. In the helium triplet state the two electrons can have the same spin 
orientation, and therefore our result represents a final rebuttal of the objection of Nesbet. 



II. CONSTRUCTION OF THE KOHN-SHAM POTENTIAL FROM ELECTRON 

DENSITY 



In the present work we apply the scheme of van Leeuwen and Baerends to construct accurate 
KS potentials for the lowest states of the helium atom. Following van Leeuwen and Baerends, we 
obtain after multiplying the KS equations 10 from the left by 4>*(v) and summing over the N 
electrons 

N 

VKs(r)p(r) = Y,[yi( r ) v2< t>iW+ s i\Mr)\ 2 ], (13) 

i=l 

where p(r) is the electron density 

N 

p(r)=]T|<Mr)| 2 . (14) 

»=i 

This leads to a self-consistency problem, which can be solved by iteration. Defining the electronic 
part of the potential, v e i(r), by 

7 

VKs{r) = -—+ v c \(r) + const., (15) 
the solution is obtained by means of the formula 

«4 fl (r) = ^«&(r) l (16) 
Po{r) 

where po{r) is the exact many-body density and Pfc(r) is the density generated with the potential 
w^j(r). This procedure is continued until certain convergence criteria are met. 



III. MANY-BODY THEORY 



The many-body electron density needed for this procedure has been evaluated by means of many- 
body perturbation technique |36| , using the nonrelativistic pair-correlation program developed by 
Salomonson and Oster [37|]. We shall briefly indicate this procedure here. 

We want to solve the Schrodinger equation 

H^ = E^ (17) 
and partition the Hamiltonian into a zeroth-order hamiltonian and a perturbation 

H = H +H'. (18) 
We start from a zeroth-order or model function ^q, which is an eigenfunction of Hq, 

H y = E y . (19) 



AE = 



FIG. 1: Upper line: Graphical representation of the pair equation (Eq. {29 )• The vertical lines represent 
the valence orbitals (double arrows) and virtual orbitals (single arrow) . The thick horizontal line represents 
^2, the dotted line the electrostatic interaction between the electrons and the box the effective two-body 
interaction W% 127H . Lower line: Graphical representation of the energy shift due to the perturbation (Eq. 
{29). 



The exact solution can be expressed 

* = Sl*o, (20) 
where 51 is the wave operator, satisfying the generalized Block equation in the linked- diagram form 

[Q,H ]P=(H'Q-nw) linked P. (21) 

Here, W is the effective interaction, in intermediate normalization (^o = P^) given by W — 
PH'QP. P is the projection operator for the model space, which in this simple case is assumed to 
contain only a single model state, ^>q. Only so-called linked diagrams will contribute according to 
the linked-diagram theorem 36] . 

Using second quantization, the wave operator can be separated into normal-ordered one-, two-,., 
body parts 

n = i + n 1 + n 2 + ... (22) 

or 

= 1+ I'M'. ,K + y {ototo,o fc } x% + ■■■ (23) 

using the sum convention, /a are the electron creation/annihilation operators, and the curly 
brackets denote the normal-ordering. The n-body part of the wave operator then satisfies the 
equation 



[n n ,H ]p = (H>n-nw) liDkc ^p 



(24) 



5 



IV. APPLICATION TO THE LOWEST STATES OF THE HELIUM ATOM 

For heliumlike systems, starting from hydrogenlike orbitals, the wave operator can be expressed 
by means of the two-body part only of the wave operator, 

Q = l + n 2 , (25) 

satisfying the 'pair equation' 

[n Jl ff ]P=(ff'fl-QW 2 ) ll)W|a P. (26) 

This equation is exhibited graphically in the upper part of Fig. ^ Here, the thick line represents 
r2 2 and the box the two-body part of the effective interaction 

W 2 = (PH'VLP) 2 . (27) 

The total energy of the system is 

E = E + AE, (28) 
where the energy shift, AE, is in this case given by 

AE = (* |-ff'n|*o> = <*o|VK 2 |*o) (29) 

and represented graphically by all 'closed' two-body diagrams, as indicated at the bottom part of 
Fig. ^ Since the final state of the ionization process is in our cases the ground state of the He 4 " 
ion, with the exact nonrelativistic energy of -2 H, the binding energy of the outermost electron 
becomes (in atomic units) 

BE = -2 - E. (30) 




RADIUS r/a 

FIG. 2: The Kohn-Sham density (dots) superimposed on the many-body density (solid line) for the helium 
ground state. 

In the present work the pair equation l|26|l has been solved, using the numerical procedure devel- 
oped by Salomonson and Oster [23, and densities for the Is 2 1 S ground state and the lowest triplet 
state, ls2s 3 S, of the helium atom have been evaluated. These densities are then used to construct 
the corresponding KS potentials, as discussed above. The wave functions obtained in this way are 
virtually exact, apart from relativistic, mass-polarization and quantum-clcctrodynamical effects. 
In a similar fashion we have also used the Hartree-Fock density to construct the corresponding KS 
potential. 
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FIG. 3: The electronic part of the Kohn-Sham potential for the helium ground state. 



It can be argued that the procedure used here corresponds to approaching the electron-density 
integral (|12fl to the electron number from below, M — ► N— ^7|> an( i hence the ionization-potential 
theorem can be tested. 

In order to achieve good accuracy, particularly for the eigenvalue of the outermost electron 
orbital, it is important to have the exact density in an accurate form and to have this density well 
reproduced by the Kohn-Sham orbitals. In the present work we have generated the many-body 
density using a large numerical grid, and the convergence criteria are set so that the Kohn-Sham 
density should not deviate from the many-body density by more than one part in 10 9 at any point. 
The convergence rate was usually quite slow, and several thousands of iterations were often needed 
to reach this level of accuracy. To improve the convergence rate and avoid 'oscillations', it was 
sometimes helpful to take some average of the last two iterations as the input for the next one. 
It is also important to keep the electronic part, v \, of the KS-potential positive at all points by 
adjusting the constant in Eq. (|15|l . After the iteration procedure was completed, the constant is 
determined so that the potential approaches zero as r — > oo. 



1.2 




FIG. 4: The Kohn-Sham density (dots) superimposed on the many-body density (solid line) for ls2s 3 S. 



As in our earlier works |37|. we apply an exponential radial grid ~ e Xi /Z, where Xi is a 
discrete linear lattice with equally spaced points with x typically ranging from a; m i n = —11 to 
about x max = 4. For the triplet state at least four different grids were used, and the results 
extrapolated in the standard way. Also the end point of the grid, x max , is varied, as it was 
found that the Kohn-Sham results were quite sensitive to that value (See Fig. UJ), most likely due 
to the fact that our pair functions are forced to be zero at the end point. This has very little 
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FIG. 5: The electronic part of the Kohn-Sham potential for the ls2s 3 S state of helium. 



effect on standard many-body calculations, but the KS eigenvalue depends strongly upon the tail 
of the density distribution and therefore more affected by the boundary condition. Hence, an 
extrapolation of x max is required in the KS case. It is likely that the results could be improved by 
more sophisticated boundary conditions. 

In the evaluation of the many-body electron density a partial- wave expansion is used |37j | , nor- 
mally up to ^ max = 10, and an extrapolation performed in the standard way. In our procedure, 
however, the Kohn-Sham potential is evaluated for successive truncations of the partial-wave ex- 
pansion, and the above-mentioned ionization-potential theorem could be tested for each truncation 
separately, as well as after the l max extrapolation. 



A. The many-body density of the helium ground state 



As a preliminary test of our procedure, we have applied this to the ground state of the helium 
atom, where the Kohn-Sham potential has previously been constructed (23L l27j . The electron 
density obtained from the KS orbitals is shown in Fig. [21 (dots) superimposed on the many-body 
density (dots). In this figure, the two densities are indistinguishable. The resulting KS potential 
is shown in Fig. [3J 

These calculations have been performed for a single grid with 201 points, with a; max = 3.6 
and ^ max = 10 without any extrapolations. The results obtained is then -0.903 7041 H for the 
KS Is eigenvalue and 0.903 7052 H for the many-body ionization energy, which verifies the above- 
mentioned ionization-potential theorem to 5-6 digits. By careful extrapolations the pair-correlation 
approach yields the value 0.903 724 39 H [37|, which agrees well with the very accurate value 
obtained by Frankowski and Pekeris [3^ and by Freund et al. [3|j of 0.903 724 377 H (uncorrected 
for relativity, mass-polarization and QED effects). 

In the corresponding calculation by Harbola [27| . an electron density taken from the litera- 
ture 40] was used. Two different configurations, ls 21 S and ls2s 1 S, respectively, were used, and 
the energy eigenvalue of the highest occupied orbital was in both cases found to be 0.899 H. 



B. The many-body density of the lowest triplet state of helium 



As mentioned, our primary goal of the present work has been to construct the Kohn-Sham 
potential from the many-body density for the lowest triplet state, ls2s 3 S, of the helium atom, 
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TABLE I: Comparison between the Kohn-Sham 2s eigenvalue and the many-body ionization energy for 
the ls2s 3 S state of helium with / max = 10 and different end points of the numerical grid. 







3-max 


KS eigenvalue 


Many-body IP 


3.8 


-0.175 228 7967 


0.175 229 3578 


4.0 


-0.175 229 1111 


0.175 229 3634 


4.2 


-0.175 229 2488 


0.175 229 3649 


4.4 


-0.175 229 3135 


0.175 229 3634 


extrapol 


-0.175 229 3630 


0.175 229 3639 



TABLE II: Comparison between the Kohn-Sham 2s eigenvalue and the many-body ionization energy for 
the ls2s 3 S' state of helium with different truncations of the partial-wave expansion. 







^max 


KS eigenvalue 


Many-body IP 


4 


-0.175 228 6206 


0.175 228 6214 


6 


-0.175 229 2341 


0.175 229 2354 


8 


-0.175 229 3366 


0.175 229 3379 


10 


-0.175 229 3630 


0.175 229 3639 


extrapol 


-0.175 229 3794 


0.175 229 3797 



and in this case we have performed extensive extrapolations, as we shall demonstrate below. 

The final KS-density for the ls2s 3 S system is shown in Fig. 0] (dots), together with the many- 
body density (solid line). The corresponding KS-potential is shown in Fig. It is interesting to 
note that the potential has a 'bump' close to the node of the outermost (valence) electron. This is 
typical of this kind of potential and is an effect of the electron self interaction (SIC) [4lL . 

1/3 

This depends approximately on p val , where p V ai is the density of the valence electron, and hence 
varies strongly near the node of the valence orbital. 

In Table [fl we show the KS 2s eigenvalue and the corresponding many-body ionization energy 
after grid extrapolation with the partial-wave expansion truncated at / max = 10 and different values 
of the grid end point, x max - This extrapolation is illustrated graphically for two values of Z max in 
Fig. It is found that the KS eigenvalue is - in contrast to the many-body ionization energy - 



x 10" 7 




FIG. 6: Zmax extrapolation of the 2s eigenvalue (rings) and the negative ionization energy (crosses) for the 
ls2s 3 5 state of helium with Z max = 5 and 10. 
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FIG. 7: i m ax extrapolation of the 2s eigenvalue and the negative ionization energy for ls2s 3 S of helium. 



quite sensitive to the end point. After the x max extrapolation, the values are found to agree to 
nine digits for each partial wave truncation. 

In order to find the "true" values of these quantities, it is necessary also to extrapolate the 
partial-wave expansion, with the result shown in Table ITU This is illustrated in Fig. These final 
values represents the nonrelativistic ionization energy, as before uncorrected for mass polarization 
as well as for relativistic and QED effects. These values agree to eight digits with the corresponding 
value 0,175 229 3782 H, obtained by Pekeris 0. 



As an additional test of the procedure described here, we have applied this also to the Hartree- 
Fock density of lowest triplet state of the helium atom. This density is generated by solving the 
standard HF equations and then inserted into the generating formula (|16fl in place of the many- 
body density. The resulting densities are quite similar to those given above, as is the resulting KS 
potential, since HF is quite a good approximation for this system. 

In the HF approximation the orbital eigenvalues correspond exactly to the corresponding ion- 
ization energies (with opposite sign), and therefore a comparison of the Hartree-Fock-Kohn-Sham 
(HFKS) 2s eigenvalue with the corresponding HF value would constitute a further test of the 
above-mentioned ionization-potential theorem. Here, we found that the agreement is extremely 
good without any x ma x extrapolation. As an illustration we give the values obtained after grid 
extrapolation for .T max = 4, where the HFKS 2s eigenvalue is -0.174 256 072 542 H and the HF 
value -0.174 256 072 544 H - an agreement to 11 digits! After complete extrapolations the result 
is -0.174 256 0724 H, which is expected to be accurate to 8-9 digits. 



We have demonstrated that it is possible to construct a local Kohn-Sham potential for the lowest 
triplet state of the helium atom with extreme accuracy. The agreement between the absolute values 
of the ionization potential and highest-lying KS orbital energy eigenvalue is verified to nine digits, 
which - as far as we know - represents by far the most accurate numerical test of the ionization- 
potential theorem performed to date. This result also clearly demonstrates that there is no conflict 



C. The Hartree-Fock density of the lowest triplet state of helium 



V. 



SUMMARY AND COMMENTS 
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between the locality theorem and the exclusion principle, as claimed by Nesbet |33|, |3J, |3 
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